---
title: "maps_iraq"
author: "MNB"
date: "1/11/2025"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



## setting libraries
```{r, include=FALSE}
library(dplyr)
library(tidyr)

library(geosphere)
library(rworldmap)
library(rgdal)
library(sf)
library(rgeos)
library(maps)
library(mapdata)
library(ggmap)
library(tmap)
library(GADMTools)
library(mapproj)
library(cartography)
library(elevatr)
library(raster)
library(sp)   
library(PBSmapping)
library(countrycode)


library(naniar)
library(zoo)
library(stringr)
library(tidytext)

```

# Remove lists 
```{r, include=FALSE}
## remove lists
rm(list=ls())
```

# Functions
```{r, include=FALSE}

'%!in%' <- function(x,y)!('%in%'(x,y))

## drop NAs
Drop_NA <- function(dataset, col_name){
  col_name <- enquo(col_name)
  dataset %>%
    drop_na((!!col_name)) %>%
    as.data.frame()
}

# drop var by name
Drop_Var <- function(dataset,col_name){
  col_name <- enquo(col_name)
    dataset %>% dplyr::select( -c((!!col_name)))
}

Nas_to_0 <- function(dataset, x) {
  
  dataset %>% mutate_at(all_of(x), ~replace(., is.na(.), 0))
  
}


```


## Import country shapes
```{r, include=FALSE}
## import shapefile country level
IRQ0_spdf <-st_read(file.path(getwd(), "gadm41_IRQ_0.shp")) 
## import shapefile admin2 level
IRQ_spdf <-st_read(file.path(getwd(),"geoBoundaries-IRQ-ADM2-all", "geoBoundaries-IRQ-ADM2.shp"))
## renemae the column for the names of districts
names(IRQ_spdf)[names(IRQ_spdf) == "shapeName"] <- "NAME_2"
## drop useless columns
IRQ_spdf<-IRQ_spdf[ , names(IRQ_spdf) %in% c("NAME_2")]

```

## setting data aid
```{r}

Iraq_aid<-read.csv(file.path(getwd(),"level_1a.csv"),)
df<-Iraq_aid
# select sectors of interest
#create short codes
df$Short_code<-c((substr(( c(df$ad_sector_codes)), 1, 2)))

## rename codes
# Define replacements using a named vector
replacements <- c(
  "10" = "Social_Infra",
  "11" = "Education",
  "12" = "Health",
  "13" = "Population_Policies",
  "14" = "Water_Supply_Sanitation",
  "15" = "Gov_Civ_Soc",
  "16" = "Soc_Infra",
  "20" = "Eco_Infra",
  "21" = "Trans_Stor",
  "22" = "Comm",
  "23" = "Energy_Gen_Supp",
  "24" = "Banking",
  "25" = "Business",
  "30" = "Production Sectors",
  "31" = "Agri_For_Fish",
  "32" = "Ind_Min_Cons",
  "33" = "Trade_policy",
  "41" = "Env_Protection",
  "42" = "Women",
  "43" = "Other",
  "51" = "Budget_Support",
  "52" = "Food_Security",
  "53" = "Commodity_Assistance",
  "60" = "Debt",
  "70" = "Humanitarian_Aid",
  "72" = "Emergency_Response",
  "73" = "Reconstruction_Relief",
  "74" = "Disaster_Prevention",
  "91" = "Administrative_costs",
  "92" = "Support_NGO",
  "93" = "Refugees_Donor",
  "99" = "Unallocated"
)

# Apply replacements using named vector and %in% operator
df$Short_code <- replacements[as.character(df$Short_code)]

# subset only projects of interest
df<-subset(df, df$Short_code=="Health"| df$Short_code=="Water_Supply_Sanitation"| df$Short_code=="Energy_Gen_Supp")

# select vars of interest
df<-df[ , !names(df) %in% c("project_location_id",
                                         "geoname_id",
                                         "place_name",
                                         "location_type_code",
                                         "gazetteer_adm_name",
                                         "is_geocoded",
                                         "project_title",
                                         "start_actual_type",
                                         "end_actual_type",
                                         "donors",
                                         "recipients",
                                         "status",
                                         "recipients_iso3",
                                         "location_class",
                                         "gazetteer_adm_code",
                                         "donors_iso3",
                                         "total_commitments",
                                         "total_disbursements")]


# ##assign 0s, NAs
List<- colnames(df)[colnames(df) %in% c("even_split_commitments", "even_split_disbursements")]
df<-Nas_to_0(df, List)

# keep only years of interest
df<-subset(df, df$transactions_start_year>2002 & df$transactions_end_year<2015)

#select projects at admin level 2 or minor
df<-subset(df, df$geographic_exactness<=3)
df<-df[!grepl("first-order", df$location_type_name),]

# drop if lat long NA
df<-Drop_NA(df, latitude)

## assign crs to coordinates
df <- st_as_sf(x = df,
               coords = c("longitude", "latitude"),
               crs = "+proj=longlat +datum=WGS84 +no_defs")

# assign to each event the polygon where it occurred
df <- st_intersection(x = IRQ_spdf, y = df)


df <- df %>% 
  group_by(NAME_2) %>%
  dplyr::summarize(
    dd = sum(even_split_disbursements, na.rm = TRUE),
    cc = sum(even_split_commitments, na.rm = TRUE)
  )

# join with full polygon set
df_aid_s = st_join(IRQ_spdf, df)
names(df_aid_s)[names(df_aid_s) == "NAME_2.y"] <- "NAME_2"
names(df_aid_s)[names(df_aid_s) == "dd"] <- "Total disbursements  (USD)"
names(df_aid_s)[names(df_aid_s) == "cc"] <- "Total commitments (USD)"

df_aid_s<-df_aid_s[ , names(df_aid_s) %!in% c("NAME_2.x")]

# ##assign 0s, NAs
List<- c("Total disbursements  (USD)", "Total commitments (USD)")
df_aid_s2<-Nas_to_0(df_aid_s, List) 

```



### Figure 2
```{r}
Colour_low<-c("#f6f6f6")
Colour_high<-c("#808080")
Colour_na<-c("#000000")

Legend_sizes=c(0, 1000000, 3000000, 6000000, 12000000, 50000000, 250000000, 500000000)

Credits<-  tm_credits(" Source: Iraq AIMS \n Geocoded Research Release, Version 1.3.1", 
             position=c("left", "bottom"),
             bg.color = "white")

Layout<-tm_layout(bg.color = "white",
                  frame = FALSE,
                  main.title= 'Iraq: aid allocation distribution',
                  main.title.position = "center",
                  main.title.size= 1.5,
                  legend.title.size = 1, 
                  legend.text.size = 0.8,
                  title.fontface="bold",
                  legend.height=10)

# Save the map as a Tiff
tiff(file.path(getwd(), "map1.tiff"),
     width = 3300,
     height = 2168,
     res = 600)

tm_shape(df_aid_s2) +
  tm_borders(lwd = 2,
             col = "#000000",
             lty = "dotted") + 
  tm_fill(col = 'Total disbursements  (USD)',
          breaks = Legend_sizes,
          palette = c(Colour_low, Colour_high)) +
  tm_shape(IRQ0_spdf) +
  tm_borders(lwd = 2,
             col = "#000000",
             lty = "solid") +
  tm_layout(legend.outside = TRUE) +
  Credits +
  Layout

dev.off()

```

# Import UCDP
```{r, include=FALSE}
UCDP<-read.csv(file.path(getwd(), "ged211.csv"), )
```
## setting data violence
```{r}

df_UCDP<-UCDP

#assign ISO3C
df_UCDP$iso3c <- countrycode(df_UCDP$country, "country.name", "iso3c")
## subset Iraq
df_UCDP<-subset(df_UCDP,df_UCDP$iso3c=="IRQ")
## subset year
df_UCDP<-subset(df_UCDP,df_UCDP$year>2002 & df_UCDP$year<2015)
#create event dummy
df_UCDP$Event_Dummy<-1
#change name conflict type
df_UCDP$type_of_violence[df_UCDP$type_of_violence=="1"] <- "state_based"
df_UCDP$type_of_violence[df_UCDP$type_of_violence=="2"] <- "non_state"
df_UCDP$type_of_violence[df_UCDP$type_of_violence=="3"] <- "one_sided"

#select state based
df_Battles_map<-subset(df_UCDP, df_UCDP$type_of_violence=="state_based")
## keep only actor b !government
df_Battles_map<-df_Battles_map[!grepl("Government", df_Battles_map$side_b),]

## keep only good precision
df_Battles_map<-subset(df_Battles_map, df_Battles_map$where_prec<=3)

# create ratio
df_Battles_map$LER<-(df_Battles_map$deaths_b+1)/(df_Battles_map$deaths_a+1)
df_Battles_map$LER<-ifelse((df_Battles_map$deaths_b+df_Battles_map$deaths_a)==0, NA,df_Battles_map$LER)

# create GTME
df_Battles_map$GTME<-(df_Battles_map$deaths_b)/(df_Battles_map$deaths_a+df_Battles_map$deaths_b)
df_Battles_map$GTME<-ifelse((df_Battles_map$deaths_b+df_Battles_map$deaths_a)==0, NA,df_Battles_map$GTME)

prj_dd<-"+proj=longlat +datum=WGS84 +no_defs"
# convert to sf object
df_Battles_map<- st_as_sf(x = df_Battles_map,
                           coords = c("longitude", "latitude"),
                           crs = prj_dd)
# join with full polygon set
df_Battles_map = st_join(IRQ_spdf, df_Battles_map)

df_Battles_map_s<- df_Battles_map %>%
                    dplyr::group_by(NAME_2) %>%
                      dplyr::summarize(`Total battles` = sum(Event_Dummy, na.rm = TRUE),
                                        LER = mean(LER, na.rm = TRUE),
                                        GTME = mean(GTME, na.rm = TRUE))



df_Battles_map_s$`Total battles`[is.na(df_Battles_map_s$`Total battles`)] = 0

```

### Figure 3
```{r}
Colour_low<-c("#f6f6f6")
Colour_high<-c("#808080")
Colour_na<-c("#000000")

Legend_sizes=c(0,1, 5, 10, 50, 100, 200, 400, 700)


Credits<-  tm_credits("Source: UCDP GED", 
             position=c("left", "bottom"),
             bg.color = "white")

Layout<-tm_layout(bg.color = "white",
                  frame = FALSE,
                  main.title= 'Iraq: battles distribution',
                  main.title.position = "center",
                  main.title.size= 1.5,
                  legend.title.size = 1, 
                  legend.text.size = 0.8,
                  title.fontface="bold",
                  legend.height=10)

# Save the map as a Tiff
tiff(file.path(getwd(), "map2.tiff"),
     width = 3300,
     height = 2168,
     res = 600)

tm_shape(df_Battles_map_s)+
  tm_borders(lwd=2,
             col = "#000000",
             lty = "dotted")+ 
  tm_fill(col='Total battles',
          breaks = Legend_sizes,
          palette=c(Colour_low,Colour_high))+
tm_shape(IRQ0_spdf)+
  tm_borders(lwd=2,
             col = "#000000",
             lty = "solid")+
  tm_layout(legend.outside = TRUE)+
  Layout+
  Credits

dev.off()


```

### Figure 4
```{r}

Colour_low<-c("#f6f6f6")
Colour_high<-c("#808080")


Legend_sizes=c(0.01, .2, .4, 0.5, .6, .7, .8, 1)

df_Battles_map_s_0s<-subset(df_Battles_map_s,df_Battles_map_s$`Total battles`==0)

## select regions
Zebra_list<-df_Battles_map_s_0s$NAME_2

# subset regions
Zebra_regions<-subset(IRQ_spdf, IRQ_spdf$NAME_2 %in% Zebra_list)

## merge polygon
polygon_sf_Zebra_regions<-(st_union(Zebra_regions$geom))

library(cartography)
# creating the zebra pattern
polygon_sf_Zebra_regions<-polygon_sf_Zebra_regions %>% hatchedLayer(mode = "sfc", pattern = "horizontal", density = 3)
polygon_sf_Zebra_regions <- st_sf(geometry = polygon_sf_Zebra_regions)


Zebra_regions_map<-tm_shape(polygon_sf_Zebra_regions) +
    tm_lines(col = "black")

df_Battles_map_s$GTME[is.na(df_Battles_map_s$GTME)] = 0

Layout<-tm_layout(bg.color = "white",
                  frame = FALSE,
                  main.title= 'Iraq: average GTME',
                  main.title.position = "center",
                  main.title.size= 1.5,
                  legend.title.size = 1, 
                  legend.text.size = 0.8,
                  title.fontface="bold",
                  legend.height=10)

# Save the map as a Tiff
tiff(file.path(getwd(), "map3.tiff"),
     width = 3300,
     height = 2168,
     res = 600)

tm_shape(df_Battles_map_s)+
  tm_borders(lwd=2,
             col = "#000000",
             lty = "dotted")+ 
  tm_fill(col='GTME',
          breaks = Legend_sizes,
          palette=c(Colour_low,Colour_high))+
tm_shape(IRQ0_spdf)+
  tm_borders(lwd=2,
             col = "#000000",
             lty = "solid")+
  tm_layout(legend.outside = TRUE)+
  Zebra_regions_map+
  tm_add_legend('line', 
                col = c("black"),
                border.col =  c("black"),
                alpha = 1,
                lwd=2,
                labels = c("NA"),
                title="",
                is.portrait = T)+
  Layout+
  Credits

dev.off()



```


## Add libraries
```{r, include=FALSE}
##
library(haven)
library(readxl)
library(foreign)
library(readr)
library(ivreg)
library(sandwich)
library(lmtest)
library(ggplot2)
# Load the required package
library(gridExtra)
# Load the required packages
library(cowplot)

library(Hmisc)
library(data.table)
library(dvmisc)
library(roll)
library(reshape)
```

## Load data
```{r}
# Import the data
data_og <- read_dta(file.path(getwd(),"Iraq_2.dta"))

```

# Define the formula for 2SLS
```{r}
formula_iv <- GTME ~ logi_SX + battles + elevation + bw_off + sh_off + border_control + mil_support + Sunni + Shia + Kurds + Turkmen + lnAid*dist_baghdad + factor(ID_side_b) |
                      logi_SX + battles + elevation + bw_off + sh_off + border_control + mil_support + Sunni + Shia + Kurds + Turkmen + factor(ID_side_b) + food_av_4_lag_6*dist_baghdad + diff_02_04*dist_baghdad


```

# Fit the 2SLS model
```{r}
data_M1 <- data_og

model_iv <- ivreg(formula_iv, data = data_M1)

# Compute clustered standard errors at the ID_name2 level
clustered_se <- vcovCL(model_iv, cluster = ~ID_name2)

# Summary with clustered standard errors
coeftest(model_iv, vcov = clustered_se)
```


## Figure 5
```{r}
data <- data_M1

dist_seq <- seq(0, 51, by = 0.01)

mean_logi_SX         <- mean(data$logi_SX, na.rm = TRUE)
mean_battles         <- mean(data$battles, na.rm = TRUE)
mean_elevation       <- mean(data$elevation, na.rm = TRUE)
mean_bw_off          <- mean(data$bw_off, na.rm = TRUE)
mean_sh_off          <- mean(data$sh_off, na.rm = TRUE)
mean_border_control  <- mean(data$border_control, na.rm = TRUE)
mean_mil_support     <- mean(data$mil_support, na.rm = TRUE)
mean_Sunni           <- mean(data$Sunni, na.rm = TRUE)
mean_Shia            <- mean(data$Shia, na.rm = TRUE)
mean_Kurds           <- mean(data$Kurds, na.rm = TRUE)
mean_Turkmen         <- mean(data$Turkmen, na.rm = TRUE)

mean_food            <- mean(data$food_av_4_lag_6, na.rm = TRUE)
mean_diff            <- mean(data$diff_02_04, na.rm = TRUE)

# Step 4: Define a reference country
reference_country <- levels(data$ID_side_b)[1]  # Use the first level as reference, or replace with a specific ID


# Step 4: Create a new data frame for prediction with CR_reign_mean varying and other variables fixed at their mean
pred_data <- data.frame(
  dist_baghdad = dist_seq,
  logi_SX = mean_logi_SX,
  battles = mean_battles,
  elevation = mean_elevation,
  bw_off = mean_bw_off,
  sh_off = mean_sh_off,
  border_control = mean_border_control,
  mil_support = mean_mil_support,
  Sunni = mean_Sunni,
  Shia = mean_Shia,
  Kurds = mean_Kurds,
  Turkmen = mean_Turkmen,
  lnAid = .01,
  food_av_4_lag_6 = mean_food,
  diff_02_04 = mean_diff,
  ID_side_b = factor(1)  # Fix ID at reference country
)

# Step 5: Predict GTME and get the standard errors for confidence intervals
predictions <- predict(model_iv, newdata = pred_data, se.fit = TRUE, vcov = clustered_se)
predictions <-as.data.frame(predictions)

# Extract the fitted values (predictions) and standard errors
pred_data$GTME <- predictions$fit
pred_data$se <- predictions$se.fit

# Step 6: Calculate the confidence intervals (e.g., 95% CI)
pred_data$ci_lower <- pred_data$GTME - 1.96 * pred_data$se
pred_data$ci_upper <- pred_data$GTME + 1.96 * pred_data$se

# Plot 1
(Plot1 <- ggplot(pred_data, aes(x = dist_baghdad, y = GTME)) +
  geom_line(color = "black") +  # Line for the predicted values
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2, fill = "#5B5B5B") +  # Shaded CI region
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  geom_rug(data = data, aes(x = dist_baghdad), sides = "b", color = "black", alpha = 0.5, inherit.aes = FALSE) +  # Rug plot on the bottom x-axis
  labs(title = "Aid = .01",
       x = "Distance",
       y = "Predicted GTME") +
  theme_light() +  # Light theme with grid lines
  theme(panel.background = element_rect(fill = "white", color = NA)) +  # Ensure white background
  scale_x_continuous(breaks = scales::pretty_breaks(5)) +  # Increase x-axis tick marks
  scale_y_continuous(limits = c(-2, 2), breaks = scales::pretty_breaks(10))    # Set y-axis limits
)

## P2, right plot
data <- data_M1

dist_seq <- seq(0, 51, by = 0.01)

mean_logi_SX         <- mean(data$logi_SX, na.rm = TRUE)
mean_battles         <- mean(data$battles, na.rm = TRUE)
mean_elevation       <- mean(data$elevation, na.rm = TRUE)
mean_bw_off          <- mean(data$bw_off, na.rm = TRUE)
mean_sh_off          <- mean(data$sh_off, na.rm = TRUE)
mean_border_control  <- mean(data$border_control, na.rm = TRUE)
mean_mil_support     <- mean(data$mil_support, na.rm = TRUE)
mean_Sunni           <- mean(data$Sunni, na.rm = TRUE)
mean_Shia            <- mean(data$Shia, na.rm = TRUE)
mean_Kurds           <- mean(data$Kurds, na.rm = TRUE)
mean_Turkmen         <- mean(data$Turkmen, na.rm = TRUE)

mean_food            <- mean(data$food_av_4_lag_6, na.rm = TRUE)
mean_diff            <- mean(data$diff_02_04, na.rm = TRUE)

# Step 4: Define a reference country
reference_country <- levels(data$ID_side_b)[1]  # Use the first level as reference, or replace with a specific ID


# Step 4: Create a new data frame for prediction with CR_reign_mean varying and other variables fixed at their mean
pred_data <- data.frame(
  dist_baghdad = dist_seq,
  logi_SX = mean_logi_SX,
  battles = mean_battles,
  elevation = mean_elevation,
  bw_off = mean_bw_off,
  sh_off = mean_sh_off,
  border_control = mean_border_control,
  mil_support = mean_mil_support,
  Sunni = mean_Sunni,
  Shia = mean_Shia,
  Kurds = mean_Kurds,
  Turkmen = mean_Turkmen,
  lnAid = 1.8,
  food_av_4_lag_6 = mean_food,
  diff_02_04 = mean_diff,
  ID_side_b = factor(1)  # Fix ID at reference country
)

# Step 5: Predict GTME and get the standard errors for confidence intervals
predictions <- predict(model_iv, newdata = pred_data, se.fit = TRUE, vcov = clustered_se)
predictions<-as.data.frame(predictions)

# Extract the fitted values (predictions) and standard errors
pred_data$GTME <- predictions$fit
pred_data$se <- predictions$se.fit

# Step 6: Calculate the confidence intervals (e.g., 95% CI)
pred_data$ci_lower <- pred_data$GTME - 1.96 * pred_data$se
pred_data$ci_upper <- pred_data$GTME + 1.96 * pred_data$se

# Plot 1
(Plot2 <- ggplot(pred_data, aes(x = dist_baghdad, y = GTME)) +
  geom_line(color = "black") +  # Line for the predicted values
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2, fill = "#5B5B5B") +  # Shaded CI region
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  geom_rug(data = data, aes(x = dist_baghdad), sides = "b", color = "black", alpha = 0.5, inherit.aes = FALSE) +  # Rug plot on the bottom x-axis
  labs(title = "Aid= 1.8",
       x = "Distance",
       y = "") +
  theme_light() +  # Light theme with grid lines
  theme(panel.background = element_rect(fill = "white", color = NA)) +  # Ensure white background
  scale_x_continuous(breaks = scales::pretty_breaks(5)) +  # Increase x-axis tick marks
  scale_y_continuous(limits = c(-2, 2), breaks = scales::pretty_breaks(10))    # Set y-axis limits
)

# Create the combined plot
combined_plot <- plot_grid(Plot1, Plot2, labels = NULL, ncol = 2)

# Create the title
title <- ggdraw() + draw_label("", fontface = 'bold', hjust = 0.5)

# Combine the title and combined plot
final_plot <- plot_grid(title, combined_plot, ncol = 1, rel_heights = c(0.1, 1))

# Add a white background to the final plot
final_plot <- ggdraw() +
  draw_grob(ggplotGrob(final_plot)) +
  theme(plot.background = element_rect(fill = "white", color = NA))

# Save the plot as a Tiff
ggsave(
  file.path(getwd(), "plot1.tiff"),
  plot = final_plot,
  width = 8,
  height = 4,
  units = "in",
  dpi = 600,
  device = "tiff"
)

print(final_plot)
```


## Figure 6
```{r}
data <- data_M1

lnAid_seq <- seq(0, 4, by = 0.02)

mean_logi_SX         <- mean(data$logi_SX, na.rm = TRUE)
mean_battles         <- mean(data$battles, na.rm = TRUE)
mean_elevation       <- mean(data$elevation, na.rm = TRUE)
mean_bw_off          <- mean(data$bw_off, na.rm = TRUE)
mean_sh_off          <- mean(data$sh_off, na.rm = TRUE)
mean_border_control  <- mean(data$border_control, na.rm = TRUE)
mean_mil_support     <- mean(data$mil_support, na.rm = TRUE)
mean_Sunni           <- mean(data$Sunni, na.rm = TRUE)
mean_Shia            <- mean(data$Shia, na.rm = TRUE)
mean_Kurds           <- mean(data$Kurds, na.rm = TRUE)
mean_Turkmen         <- mean(data$Turkmen, na.rm = TRUE)

mean_food            <- mean(data$food_av_4_lag_6, na.rm = TRUE)
mean_diff            <- mean(data$diff_02_04, na.rm = TRUE)

# Step 4: Define a reference country
reference_country <- levels(data$ID_side_b)[1]  # Use the first level as reference, or replace with a specific ID


# Step 4: Create a new data frame for prediction with CR_reign_mean varying and other variables fixed at their mean
pred_data <- data.frame(
  lnAid = lnAid_seq,
  logi_SX = mean_logi_SX,
  battles = mean_battles,
  elevation = mean_elevation,
  bw_off = mean_bw_off,
  sh_off = mean_sh_off,
  border_control = mean_border_control,
  mil_support = mean_mil_support,
  Sunni = mean_Sunni,
  Shia = mean_Shia,
  Kurds = mean_Kurds,
  Turkmen = mean_Turkmen,
  dist_baghdad = 35,
  food_av_4_lag_6 = mean_food,
  diff_02_04 = mean_diff,
  ID_side_b = factor(1)  # Fix ID at reference country
)

# Step 5: Predict GTME and get the standard errors for confidence intervals
predictions <- predict(model_iv, newdata = pred_data, se.fit = TRUE, vcov = clustered_se)
predictions<-as.data.frame(predictions)

# Extract the fitted values (predictions) and standard errors
pred_data$GTME <- predictions$fit
pred_data$se <- predictions$se.fit

# Step 6: Calculate the confidence intervals (e.g., 95% CI)
pred_data$ci_lower <- pred_data$GTME - 1.96 * pred_data$se
pred_data$ci_upper <- pred_data$GTME + 1.96 * pred_data$se

# Plot 1
(Plot1 <- ggplot(pred_data, aes(x = lnAid, y = GTME)) +
  geom_line(color = "black") +  # Line for the predicted values
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2, fill = "#5B5B5B") +  # Shaded CI region
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  geom_rug(data = data, aes(x = lnAid), sides = "b", color = "black", alpha = 0.5, inherit.aes = FALSE) +  # Rug plot on the bottom x-axis
  labs(title = "Distance= 35",
       x = "Aid",
       y = "") +
  theme_light() +  # Light theme with grid lines
  theme(panel.background = element_rect(fill = "white", color = NA)) +  # Ensure white background
  scale_x_continuous(breaks = scales::pretty_breaks(5)) +  # Increase x-axis tick marks
  scale_y_continuous(limits = c(-2, 2), breaks = scales::pretty_breaks(10))    # Set y-axis limits
)


## P1, left plot
data <- data_M1

lnAid_seq <- seq(0, 4, by = 0.02)

mean_logi_SX         <- mean(data$logi_SX, na.rm = TRUE)
mean_battles         <- mean(data$battles, na.rm = TRUE)
mean_elevation       <- mean(data$elevation, na.rm = TRUE)
mean_bw_off          <- mean(data$bw_off, na.rm = TRUE)
mean_sh_off          <- mean(data$sh_off, na.rm = TRUE)
mean_border_control  <- mean(data$border_control, na.rm = TRUE)
mean_mil_support     <- mean(data$mil_support, na.rm = TRUE)
mean_Sunni           <- mean(data$Sunni, na.rm = TRUE)
mean_Shia            <- mean(data$Shia, na.rm = TRUE)
mean_Kurds           <- mean(data$Kurds, na.rm = TRUE)
mean_Turkmen         <- mean(data$Turkmen, na.rm = TRUE)

mean_food            <- mean(data$food_av_4_lag_6, na.rm = TRUE)
mean_diff            <- mean(data$diff_02_04, na.rm = TRUE)

# Step 4: Define a reference country
reference_country <- levels(data$ID_side_b)[1]  # Use the first level as reference, or replace with a specific ID


# Step 4: Create a new data frame for prediction with CR_reign_mean varying and other variables fixed at their mean
pred_data <- data.frame(
  lnAid = lnAid_seq,
  logi_SX = mean_logi_SX,
  battles = mean_battles,
  elevation = mean_elevation,
  bw_off = mean_bw_off,
  sh_off = mean_sh_off,
  border_control = mean_border_control,
  mil_support = mean_mil_support,
  Sunni = mean_Sunni,
  Shia = mean_Shia,
  Kurds = mean_Kurds,
  Turkmen = mean_Turkmen,
  dist_baghdad = 5,
  food_av_4_lag_6 = mean_food,
  diff_02_04 = mean_diff,
  ID_side_b = factor(1)  # Fix ID at reference country
)

# Step 5: Predict GTME and get the standard errors for confidence intervals
predictions <- predict(model_iv, newdata = pred_data, se.fit = TRUE, vcov = clustered_se)
predictions<-as.data.frame(predictions)

# Extract the fitted values (predictions) and standard errors
pred_data$GTME <- predictions$fit
pred_data$se <- predictions$se.fit

# Step 6: Calculate the confidence intervals (e.g., 95% CI)
pred_data$ci_lower <- pred_data$GTME - 1.96 * pred_data$se
pred_data$ci_upper <- pred_data$GTME + 1.96 * pred_data$se

# Plot 1
(Plot2 <- ggplot(pred_data, aes(x = lnAid, y = GTME)) +
  geom_line(color = "black") +  # Line for the predicted values
  geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper), alpha = 0.2, fill = "#5B5B5B") +  # Shaded CI region
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  geom_rug(data = data, aes(x = lnAid), sides = "b", color = "black", alpha = 0.5, inherit.aes = FALSE) +  # Rug plot on the bottom x-axis
  labs(title = "Distance= 5",
       x = "Aid",
       y = "Predicted GTME") +
  theme_light() +  # Light theme with grid lines
  theme(panel.background = element_rect(fill = "white", color = NA)) +  # Ensure white background
  scale_x_continuous(breaks = scales::pretty_breaks(5)) +  # Increase x-axis tick marks
  scale_y_continuous(limits = c(-2, 2), breaks = scales::pretty_breaks(10))    # Set y-axis limits
)



# Create the combined plot
combined_plot <- plot_grid(Plot2, Plot1, labels = NULL, ncol = 2)

# Create the title
title <- ggdraw() + draw_label("", fontface = 'bold', hjust = 0.5)

# Combine the title and combined plot
final_plot <- plot_grid(title, combined_plot, ncol = 1, rel_heights = c(0.1, 1))

# Add a white background to the final plot
final_plot <- ggdraw() +
  draw_grob(ggplotGrob(final_plot)) +
  theme(plot.background = element_rect(fill = "white", color = NA))

# Save the plot as a Tiff
ggsave(
  file.path(getwd(), "plot2.tiff"),
  plot = final_plot,
  width = 8,
  height = 4,
  units = "in",
  dpi = 600,
  device = "tiff"
)


print(final_plot)
```



```{r}

```


